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Abstract 



We develop a numerical code to compute gravitational waves induced by a particle 
moving on eccentric inclined orbits around a Kerr black hole. For such systems, the 
black hole perturbation method is applicable. The gravitational waves can be evaluated 
by solving the Teukolsky equation with a point like source term, which is computed from 
the stress-energy tensor of a test particle moving on generic bound geodesic orbits. In 
our previous papers, we computed the homogeneous solutions of the Teukolsky equation 
using a formalism developed by Mano, Suzuki and Takasugi and showed that we could 
compute gravitational waves efficiently and very accurately in the case of circular orbits 
on the equatorial plane. Here, we apply this method to eccentric inclined orbits. The 
geodesies around a Kerr black hole have three constants of motion: energy, angular 
momentum and the Carter constant. We compute the rates of change of the Carter 
constant as well as those of energy and angular momentum. This is the first time 
that the rate of change of the Carter constant has been evaluated accurately. We also 
treat the case of highly eccentric orbits with e = 0.9. To confirm the accuracy of our 
codes, several tests are performed. We find that the accuracy is only limited by the 
truncation of £-, k- and n-modes, where I is the index of the spin-weighted spheroidal 
harmonics, and n and k are the harmonics of the radial and polar motion, respectively. 
When we set the maximum of i to 20, we obtain a relative accuracy of 10~ 5 even in 
the highly eccentric case of e = 0.9. The accuracy is better for lower eccentricity. 
Our numerical code is expected to be useful for computing templates of the extreme 
mass ratio inspirals, which is one of the main targets of the Laser Interferometer Space 
Antenna (LISA). 
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§1. Introduction 



Gravitational radiation is one of the most important predictions resulting from general 
relativity. The presence of gravitational radiation has been proved indirectly through its 
effect on the orbital period of the Hulse- Taylor binary pulsar. ^ Owing to advances in modern 
technology, the era of gravitational wave astronomy has almost arrived. Several ground-based 
interferometric gravitational wave detectors have been operated 2 ** -5 ** for several years. Next- 
generation ground-based interferometers such as advanced LIGO, advanced VIRGO and 
LCGT, are planned, and will be started in the near future. Research and development studies 
of a space-based gravitational wave observatory project, the Laser Interferometer Space 
Antenna (LISA), 6 ) are rapidly progressing. There are also proposals for laser interferometer 
gravitational wave antennas in space such like a DECihertz Interferometer Gravitational 
wave Observatory (DECIGO) 7 * 1 and a Big Bang Observer (BBO). 8 * 1 Those will be sensitive 
to frequencies of 10 -2 < / < 1 Hz. 

One of the most promising sources of gravitational waves that can be detected by LISA 
is a compact star orbiting a supermassive black hole. Observing gravitational waves from 
this type of binary system, i.e., an extreme mass ratio inspiral (EMRI), we may be able to 
obtain information on the central black hole's spacetime such as the mass, spin of the black 
hole and the mass distribution of compact objects in the center of the galaxy. To extract 
the physical information of EMRI from observational data obtained by the detectors, we 
have to compute the theoretical waveforms with a phase accuracy within one cycle over the 
total number of cycles. LISA is sensitive to gravitational waves around 10 -2 Hz. When the 
observation of LISA is performed for one year, the total number of cycles of waves is ~ 10 5 . 
Thus, to analyze one year of data obtained from LISA, we need theoretical waveforms that 
are accurate to 10 -5 . 

The dynamics of EMRI is accurately modeled as a point particle of small mass moving 
around a Kerr black hole. Therefore, gravitational waves from EMRI can be evaluated using 
black hole perturbation theory, which was originally developed as a metric perturbation 
theory for a black hole spacetime. For nonrotating (Schwarzschild) black holes, a single 
master equation for the metric perturbation was derived by Regge and Wheeler for the 
odd-parity parts, 9 '' and later by Zerilli for the even-parity parts. 10 * 1 These equations are 
remarkably simple, separable, hyperbolic equations with a potential term. For rotating black 
holes, at present, there are no simple decoupled equations for metric coefficients. Instead 
the perturbed geometry must be analyzed by the equations derived by Teukolsky using 
gauge-invariant variables corresponding to some tetrad components of the perturbed Weyl 
curvature. n ) 
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Using this Teukolsky formalism, there have been many numerical computations of gravi- 
tational waves induced by a point particle. See Chandrasekhar 12 ) and Nakamura et al., 13 ^ for 
reviews and for references on earlier papers. For simple orbits such as circular or equatorial 
orbits around a black hole, an accuracy of 10~ 5 has been achieved, which may be sufficient 
to detect gravitational waves. Computations of gravitational waves induced by a particle 
moving on eccentric nonequatorial orbits are now available. 14 '' However, a large amount of 
computation time is needed and the results have only a few orders of accuracy. 

Furthermore, it has been pointed out that EMRIs radiating gravitational waves within 
the LISA band have high eccentricities. 15 '' So far, the gravitational energy flux has been 
computed only for orbits with eccentricity less than 0.7. This is primarily due to the low nu- 
merical accuracy and high computational time at the high-frequency modes. In the Teukolsky 
formalism in the frequency domain, when we treat a large eccentricity and large inclination 
angle, we need to compute a large number of harmonics corresponding to the radial and 
polar motion. Thus, the computation time becomes a serious problem. It is thus desirable 
to develop more efficient and more accurate numerical codes. 

In the Kerr spacetime, there are three constants of motion of a geodesic: the energy, 
angular momentum and Carter constant. Using the conservation law, we can evaluate the 
rates of change of the energy and angular momentum of a particle due to the emission of 
gravitational waves. In contrast, the rate of change of the Carter constant cannot be derived 
from the conservation law. Mino proposed a method to evaluate the average rates of change 
of the three constants including the Carter constant 16 ^ under an adiabatic approximation. 
He showed that the average rates of change can be evaluated using the radiative field instead 
of the retarded field. Using Mino's method, Drasco et a/. 17 -* derived simplified version of his 
formula for the case of a scalar field. Sago et a/. 18 -* developed Mino's method further, giving 
a simplified version of his formula. Applying their new scheme, they gave explicit analytic 
formulas for the rates of change of constants when a particle moves on slightly eccentric and 
slightly inclined orbits. 19 -* Ganz et al. extended this computation to the case when a particle 
moves on slightly eccentric and arbitrarily inclined orbits. 20 ^ However, these two results are 
based on the assumption of small eccentricity and the post-Newtonian approximation to the 
order 0((v/c) 5 ). 

In this paper, we compute the rates of change of the three constants of motion, including 
the Carter constant, induced by a particle moving on eccentric and inclined orbits around a 
Kerr black hole without assuming small eccentricities, small inclination angles or low velocity. 
This is the first time that the rate of change of the Carter constant has been computed 
accurately using an adiabatic approximation. Drasco and Hughes 21 ^ computed the rate of 
change of the Carter constant assuming that the inclination angle does not change. However, 
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this assumption holds only approximately, rather than exactly. 

To treat cases of large eccentricity and inclination angle, we introduce various methods. 
The numerical method used to compute the homogeneous solution of the Teukolsky equation 
is based on that of Fujita and Tagoshi, 22 )' 23 ) in which the Mano-Suzuki-Takasugi formalism 24 ) 
is used. However, since only circular, equatorial orbits are treated in these works, we extend 
the method so that we can compute the homogeneous solutions efficiently in eccentric and 
inclined cases. For the source term, we introduce analytical expressions for the radial and 
polar motion to achieve high accuracy. 

The paper is organized as follows. In § [21 we summarize the details of the method 
for computing the gravitational waves from EMRIs using black hole perturbation theory. 
In § El we explain our numerical methods. In § HJ we discuss the peaks of radial and 
polar modes briefly then we derive the rates of change of the three constants of motion 
for highly eccentric orbits. Next, we verify our code using the Schwarzschild case and by 
comparison with analytical post-Newtonian expressions. Section is devoted to a summary 
and discussion. In the Appendices, some detailed formulas are given. Throughout this paper 
we use units with c = G = 1. 

§2. Formulas for the rate of change due to gravitational wave emission 

In the Teukolsky formalism, the gravitational perturbation of a Kerr black hole is de- 
scribed in terms of the Newman-Penrose variables, $o an d #4, which satisfy the master 
equation. The Weyl scalar 1^4 is related to the amplitude of the gravitational wave at infin- 
ity as 



The master equation for ^4 can be separated into radial and angular parts if we expand $4 
in Fourier harmonic modes as 



where p = (r — mcos0) _1 , the angular function _ 2 S^(9) is the spin-weighted spheroidal 
harmonic with spin s = —2, and M and aM are the mass and angular momentum of the 
black hole, respectively. The radial function Re m u{r) satisfies the radial Teukolsky equation, 



#4 — > -(h+ — ih x ), for r — > 00. 



(2d) 




(2-2) 




(2-3) 



where A = r 2 — 2Mr + a 



2 . The potential term V(r) is given as 
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where K = (r 2 + a 2 )u — ma and A is the eigenvalue of ^S^iO). 

We solve the radial Teukolsky equation by using the Green function method. The solution 
of the Teukolsky equation with a purely outgoing property at infinity and a purely ingoing 
property at the horizon becomes 

R^(r) = ^— Uz.(r) [ dr'Rl w T lmul A~ 2 

VVlmuj I J r+ 

+*ft» [ dr'RZJ^A- 2 } , (2-5) 
where the Wronskian Wi muJ is given as 

W lmui = 2zooCfZ s BZ^ (2-6) 

and where -R^" P ( r ) satisfy ingoing/outgoing wave conditions at the horizon/infinity. The 
asymptotic forms of Rfy^(r) are expressed as 



DID 
1X1 Imw 



Dtrans A 2 „-iPr* f„„ „ . „ 

a imw n e ior r ^ r+, 

r 3 pref „iu>r* I r — 1 Dine p -iur* f nr „ . 



^UP iPr* I A2^fref p -iPr* f 

Anw » r 3 ^e ,ur * forr^oo, 



where P = u — ma/2Mr + and r* is the tortoise coordinate defined as 

2Mr + , r — r i 2Mr_ , r — r_ . 
r*=r + ^ln ± In , 2-8 



with r± = M ± VM 2 - a 2 . 

The asymptotic property of the solution at the horizon is expressed as 

Dtrans /\ 2 „—iPr* poo 

Re m .(r - r + ) = ^ trans / dr^T^" 3 ee ^V^*. (2-9) 
The solution at infinity is expressed as 

RiUr - oo) = — — - / dr'Rt^T^A- 2 = Z^rV^. (2-10) 
Using the formula of the source term T^ ma; , 25 ^ Z^^ are expressed as 

Dtrans roc 



o • ,/^trans Dine 

^ = f° ^-^iLKi), *(*)], (2-11) 
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where 



^tmui — [Rlmuj {AnnO + AfnnO + Anmo} 

dRT., r . . , d 2 R 



<™ f A _ 1 j Irrea; 4 _ _ 

1 n mnl "T ^mml J T ^2^. ^ 1 -mr?i2 



^mw ~~ [^fmu {^nn.O + A mn Q + v4 mm o} 



J pin j2 pin 

url lmLu (a \ A _ _ \ 1 ^mi 4 _ 

l^mnl ~r ^mmi / ~r ,9 J iram2 



dr err 



r=r(t),0=0(t) 



r=r(t),8=8(t) 



(2-12) 



where A nn0 and other terms are given in Appendix lAl 

The function 2^^[r(t), #(£)] is constructed from the source term of the Teukolsky equa- 
tion and depends on the orbital worldline of the star perturbing the black hole spacetime. If 
the trajectory of a compact star is eccentric and inclined from the equatorial plane, it is not 
easy to evaluate Eq. (12- 111) . This is because the radial and polar motion are coupled in the 
observer time variable t. This problem can be solved by introducing a new time variable A 
defined as dX = dr/S, where E = r 2 + a 2 cos 2 #. 16 )' 26 ) The geodesic equations become 

(S) = [( r2 + a2 ) S - aC z} 2 - A i r2 + ( C z- aS ) 2 + C ]= R ( r )> (2' 13 ) 



dcos9^ " 



dX 



where 



= C-(C + a 2 (l - £ 2 ) + C 2 Z ) cos 2 9 + a 2 (l - £ 2 ) cos 4 9 

= 0(cos9), (2-14) 

^ = $ I (r)+$ g (cos9)-a£, (2-15) 
dX 

^L = T r (r)+T e (cos9)+aC z , (2-16) 
dX 



<2> r (r) = - [£(r 2 + a 2 ) - a£,] , ^(cos „ 

ZJ L J 1 — cos z 9 

2 2 

T r (r) = T —^- [£(r 2 + a 2 ) - aC z ] , T e (cos9) = -a 2 £(l - cos 2 #), 



A 

and £, L z and C are the energy, the z-component of the angular momentum and the Carter 
constant per unit mass, respectively. The equations of radial and polar motion are decoupled 
in Eqs. ( I2-13P and (I2-I4p . For the bound orbits, r(A) and 9(X) become periodic functions that 
are independent of each other. The fundamental periods for the radial and polar motion, A r 
and Ag, are defined as 

n r max dr t A r C086 ™» dcos9 , n1f . 

A r = 2 __ , A = A , 2-17 
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The angular frequencies of the radial and polar motion become 

~ 27r ~ 27r , « 

r, = -, r 9 = -. (2-18) 

Explicit expressions for T r and Tq are given in Appendix [Bj 

We define the angle variables as w r = T r X and wg = TgX. The functions that depend 
only on r or 9 become periodic functions with respect to w r or wg, respectively, with period 
2tt. 

We expand the right-hand sides of Eqs. (12- 151) and (12- 161) into Fourier series, 

^ = ^V- ar ^ A , (249) 

k,n 

^ = V^ n e-^ A e-^ A , (2-20) 
dX J 

k,n 

where 

T k , n = -^— 2 ^ dw r dw e (T r (r)+T e (cos9) + aC z )e lkWr e mw \ (2-21) 
{ Z7T ) Jo Jo 

^M = 7Tl2 /"^ [" dw ($ r (r)+$e(cose)-aS)e ikw re inw *. (2-22) 
l^J Jo JO 

(2-23) 

Since Tk } n = and <Pk,n = in the case of k ^ and n ^ 0, we have 



dt 
dX 



dX 



where 



r + T k,oe~ ikWr + T o,ne~ mWe , (2-24) 

fc^O n^O 

r = T 00 

= T tM + r tW + aC z , (2-25) 

^ + J2®k,oe- ikvh + J2®o,ne- inWe , (2-26) 

fc^O n.^0 

Y<t> = ^oo 

= V)+T^)-a£, (2-27) 



2 /»27r /»27r 

^tM = 7T / duyT n = 7T dw e T e , 

27T ,/ 27r y 

r 0M = T" / ?>) = 77- / rfwe^e- (2-28) 
^ Jo ^ Jo 



8 



Since Tk, n = for k ^ and n ^ 0, we need to compute T^o, T 0>n , @k,o and @o,n, which can 
be easily obtained from a one-dimensional Fourier transformation. We obtain the functions 
t(X) and <f)(X) from the following formulas: 



t(\) = rx + J2 



k^O 



kT r nTa 



-invjQ 



nTa 



(2-29) 
(2-30) 

The two variables, J 1 and 7^, represent the average rates of change of t and as functions 
of A, respectively. 

When we consider the bound orbits of a point particle using A, the amplitude of the 
partial wave Z^^, defined in Eq. (12-lip . can be expanded by using a Fourier series as 



(2-31) 



where 



700,H 
' imkn 



(27T) 



2tt 



2- 



and 



^ m fcn = (m^ + kT e + nT r )/r . 
Using these functions, the gravitational waveform at infinity is expressed as 



|E, 

imkn rnkn 



imkn —2 £m V / r fl^mkn{ T * ~ t)+im<j> 



(2-32) 



(2-33) 



(2-34) 



Moreover, the time-averaged rates of change for the three constants of motion due to the 
emission of gravitational waves are expressed as 16 -*' 18 -'' 19 '' 

d£' 
~dt 

d£ 2 



H 1 = -^ 



imkn mkn 



yoo 
^imkn 



^~ &£mkn 



7 K 

imkn 



dt 

dC 
~dt 
dQ 
dt 



m 



imkn mkn 



yoo 
imkn 



&lmkn 



m 

' Imkn 



dQ 
dt 

2T t (r) 



2(a£ - C z )(a 



d£_ 
~dt 



dC z 
dt 



(2-35) 
(2-36) 
(2-37) 



dS 
~dt 



dCz 
dt 



_ 9 V TilC r 

imkn mkn 



yoo 
imkn 



^imkn 



7 H 

imkn 



(2-38) 
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Fig. 1. Plots of generic bound geodesic orbits with eccentricity e = 0.7, semilatus rectum p = 10M 
and inclination angle mc = 45°. The black hole's spin is set to a = 0.9A/. The left figure 
is expressed using a Cartesian coordinate system. The right figure is expressed by using a 
corotational system, defined as Eq. ()3-2j) . 
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Fig. 2. The left figure is the projected image of right figure in Fig. Q3 The right figure is the 
projected image of the same figure to the Z = Y tan#i nc plane. Here the coordinates (X' ,Y') 
are defined as X' = X and Y' = Y cos $inc sin $inc ■ The left figure shows that the orbit 
almost stays on the plane inclined by 6 mc from the equatorial plane. 



where 



256(2Mr + ) 5 P(P 2 + 4e 2 )(P 2 + 16e 2 )u;^ n 
rvrs ■ 



(2-39) 



where e = \J M 2 - a 2 /AMr + , P = u mkn -ma/(2Mr + ), r + = M+^/M 2 - a 2 is the location of 
the event horizon and Cj^ kn is the Teukolsky-Starobinsky constant. 27 ' Here (■ ■ • ) represents 
the time average. 
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§3. Numerical methods 



In this section, we discuss the numerical method used to compute Z^[ kn , Eq. (12-321) . The 
numerical computation is divided into four parts: geodesies, homogeneous solutions, the 
integration of Eq. (!2-32[) and the mode summation in Eqs. (12 351) (12 381) . In the following 
subsections, we explain the method for obtaining the geodesic and homogeneous solutions. 
The numerical integration of Eq. (12-321) and the mode summation are explained in § HI 

3.1. Geodesies 

Geodesies around a Kerr black hole are completely specified by the three constants of 
motion, (£, C z , C). Instead of using these three constants, it is convenient to introduce three 
orbital parameters: eccentricity e, semilatus rectum p and inclination angle 6*i nc which are 
defined in Eq. ( 1B-3I) . There is a one-to-one correspondence between the orbital parameters 
(p,e,9 inc ) and (S,C Z , C) 21 ^ 

Our method of computing geodesic motion is summarized as follows. We specify three 
orbital parameters, (p, e, 9 inc ). We first compute the constants of motion, £, C z and C, which 
correspond to (p, e, 9 inc ). We then compute the orbital frequencies of the radial and polar 
motion, T r and Tg. We find that the solutions of the r and 9 components of the geodesic 
equations are expressed in terms of Jacobi elliptic functions, which can be easily computed 
numerically (see Appendix [B]) . This enables us to evaluate the bound orbits accurately. We 
have not confirmed whether this method is faster than the numerical integration methods 
when the same accuracy is required. However, since we need to compute the geodesic 
motion only once, its computation time is negligible compared with the total computation 
time. Thus, we adopt this method. The r and 9 components of coordinates and velocity of 
the particle, r(w r ), [dr / d\](w r ) , [cos 9} (wg) and [d cos 9 /d\](wg), are respectively expressed 
using w r and wg analytically. 

The t and components of the solution of the geodesic equations can be obtained using 
the Fourier series expansion, Eqs. (I2-29P and f)2-30p . We truncate the Fourier series expansion 
at a certain frequency. Since the convergence of the Fourier series is rapid, as mentioned in 
Ref. 21), this does not cause serious problems. Explicit expressions for orbital frequencies 
and the solutions of radial and polar motion are given in Appendix [B] 

Here we demonstrate a geodesic orbit computed using our code. We start an orbit at 
(t, r, 9, 0) = (0, p/(l + e), 0, 0), and end it after several oscillation periods of radial and 
polar motion. In Fig. [Tj we plot the orbit on two different coordinate systems. The left 
figure of Fig. [T]is expressed using the Cartesian coordinate system (t, x, y, z) defined as 

t = t, x = r sin 6* cos0, y = r sin 6* sin 0, z = rcos9. (3-1) 
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This figure indicates the well-known fact that generic geodesic orbits around a Kerr black 
hole are complicated. To understand the physical meaning of the orbital parameters p, e 
and 8- mc , it is useful to express the orbit using the corotational coordinate system (t, X, Y, Z) 
defined as 

t = t, X = xcos(T prc A) - ysin(T prc A), Y = xsin(T pre A) + y cos(T prc A), Z = z, (3-2) 

where the precession frequency, T pre , is defined as T pre := Tq — T^. The right figure in Fig. CD 
is expressed in the corotational coordinate system. The left figure in Fig. [2] is a plot of the 
orbits projected onto the X = plane. This figure shows that the particle almost remains 
on the Z = Y tan# inc plane. In fact, the orbit around a Schwarzschild black hole is exactly 
on this plane. On the other hand, when the parameter p becomes small and the particle 
approaches the horizon, it becomes impossible to define such an approximate orbital plane. 
In such a case, it is not appropriate to call the parameter 9\ nc the inclination angle. 

The right figure in Fig. [2] is the projected image onto the Z = Y tan 9 inc plane. The 
minimum and maximum distances from the origin, r min and r max , are r min = p/(l + e) and 
r max = p/(l — e). In this sense, e and p are called the eccentricity and semilatus rectum, 
respectively. 

3.2. Homogeneous solutions 

Mano, Suzuki and Takasugi 24 -* (MST) formulated a method to express homogeneous so- 
lutions of the Teukolsky equation, R^^ and R^,, in series of hypergeometric functions or 
Coulomb wave functions. Fujita and Tagoshi 22 -'' 23 '' applied the MST method to the numeri- 
cal computation of gravitational waves from a particle moving on circular, equatorial orbits 
around a black hole. To read off the asymptotic amplitudes such as B mc and £> rcf , we usually 
have to numerically integrate the homogeneous Teukolsky equation from r + to a large radius. 
In the MST formalism, on the other hand, the analytical forms of these asymptotic ampli- 
tudes of the homogeneous Teukolsky equation are given. We can evaluate the asymptotic 
amplitude very accurately. Furthermore, since the convergence of the series of hypergeomet- 
ric functions is fast, we can obtain the homogeneous solutions themselves very accurately. In 
this paper, since we treat eccentric orbits, we need to evaluate the homogeneous solutions at 
many radial points. Although the convergence of the series of hypergeometric functions or 
Coulomb wave functions is very fast, the computation time required to evaluate these series 
is not negligible. We thus use the method of successive Taylor series expansions. 

We first compute a homogeneous solution at a radius r, where r min < r < r max , using the 
series of hypergeometric functions. This gives a very accurate boundary value of the solution. 
We then compute the homogeneous solution at r + h using the Taylor series expansion 
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around r. The nth derivative of the homogeneous solution at radius r is evaluated from the 
recurrence relations derived by differentiating the homogeneous Teukolsky equation. The 
relative error of the Taylor series expansion is estimated as (d n \r) / 'dr n )h n j 'n\ / 'RfJ^ \r) . 
We adjust n and h to obtain a relative error ~ 10 -15 . Next, we repeat the Taylor series 
expansion around r + h to obtain the homogeneous solution near r + h. In this way, we 
obtain -R^" P (r) for all radii. 

The Taylor series expansion is used to solve ordinary differential equations when high 
accuracy is required. The Taylor series method in Ref. 29) is faster than the standard 
numerical integration methods when the required accuracy is better than 10~ 6 . We adopt 
the Taylor series method since we aim to develop a highly accurate code. We confirmed that 
this method gives accurate results and is much faster than using the hypergeometric series 
expansion. However, we have not compared the computing time of the Taylor series method 
with that of the standard numerical integration method by setting the same accuracy. The 
use of the numerical integration method together with the MST formalism may be useful for 
reducing the computation time when the required accuracy is not very stringent. 

The spin- weighted spheroidal harmonics, _25^(6 I ), are evaluated using a series of Jacobi 
polynomials. The details of the numerical method are described in Ref. 22). Although it 
should be possible to use the Taylor series method for the spin-weighted spheroidal harmon- 
ics, we have not attempted this yet. 



§4. Results 



4.1. Energy spectrum 

In this section, we compute the rates of change of the three constants of motion, including 
the Carter constant. 

We use the trapezium rule to compute the integral, Eq. (12-321) . This is because the 
integrand is a periodic function. It is well known that if a periodic function is integrated 
over one period, the trapezium rule is a suitable method for numerical integration. 

The MST formalism can be applied only for the case that the frequency is positive, u > 0. 
The mode summation in Eq. (I2-35H is performed using the following equations: 



OO OO lie \ OO 



j£ V OO oo I OO OO I 

f) = 2 E E E E(^ ) • <«> 

al I GW p = n m= _p fc= _™ „.=„.„ \ al I GW 



=2 m=— £ k=—oo n=no 
oo I oo oo/,.-, \H 



\ H oo £ oo oo / , p 

|) = 2 E E E E(%) • («) 

ah I GW f=2 m= -e fc=-oo n=n X " ' GW 
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(4-3) 
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mkn 
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imkn 
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(4-4) 



where no is the minimum integer such that mT^, + /cT# + rioT r > holds for each m and k. 
We also define the intermediate modal energy flux for later use: 



The behavior of the spectrum in terms of n is complicated. In Fig. [3], we show the modal 
energy flux at infinity, ^d£i m ^ n j dt)^ w , as functions of n in the Schwarzschild case and in the 
cases of e = 0.1, 0.5, 0.7 and 0.9. Figure H]is the same figure in the Kerr case with a = 0.9M. 
From these figures, we find that the number of peaks of (d£i m k n / dt)^ w is roughly £, and 
that the value of n at the highest peak of (dSt m kn/ <^)gw becomes larger as either £ or the 
eccentricity e becomes larger. For example, the location of the peak of {d£i m k n / ^)gw * s 
approximately n = 700 when e = 0.9, £ = m = 20 and k = 0. The shapes of the spectra in 
the Schwarzschild and Kerr cases are qualitatively very similar. 

In practical computations, we have to truncate the mode summation. We set the target 
accuracy of the mode summation with respect to n and k to be 10~ 10 . It is useful if we 
know the location of the highest peak of (d£i m kn/ ^)gw before the computation to reduce 
the computational time. However, it is difficult to derive an analytical expression for such 
a location. In this work, we adopt the following procedure to determine the range of the 
n-mode summation. First, we compute (d£i m k n / dt)cw ^ or ^ = 171 = ^ anc ^ ^ = for a 
wide range of n starting from n = uq to obtain sufficient coverage of the n-mode. In this 
computation, we obtain the location of the peak n = n p for i = m = 2 and k = 0. For 
other £ = 2 modes, since it is expected that the peak value of n is not significantly different 
from that for £ = m = 2 and k = 0, we first search for the peak near n p . We then compute 
(dSimkn/ ^)gw starting from the new peak of n, and sum {d£i mkn / dt)^ w until the total 
flux converges with an accuracy of 10~ 10 . For modes (£, m, k) with £ > 2, we basically repeat 
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(4-6) 




(4-7) 



(4-8) 
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the above procedure. We search for the peak near the peak for (£, m,k) = (£ — 1, £ — 1, 0). 
Starting from this location of n, we compute and sum {dSe m k n / ^)q W until the total flux 
converges. 

We show the spectrum of {d£i mk / dt)^ w as functions of k in Figs. \5\ and M in the 
Schwarzschild case, and in Figs. [7] and [8] in the Kerr case with a = 0.9M. Figures [5] and [7] 
show the spectrum for the inclination angle #i nc = 20°, and Figs. [6] and [8] show the spectrum 
for the high inclination angles #; nc = 70° and 80°, respectively. We find that the peak of the 
fc-mode usually exists at approximately £ — m except for the m < modes in the case of a 
low inclination angle. We also find that when the inclination angle is 20°, the peak value of 
the modal energy flux becomes smaller when m changes from £ to —£. On the other hand, 
when the inclination angle is large, i.e., 9 mc = 70° or 80°, the peak value is largest when 
m — 0. Similarly to the ra-mode, the difference between the Schwarzschild and Kerr cases is 
not very large, but the spectra of the Kerr cases are broader than those of the Schwarzschild 
cases. By taking this behavior into account, it is straightforward to truncate the k-mode 
summation, which guarantees the accuracy of the total flux of 10~ 10 . 

In Figs. M and [10], we show the modal energy flux at the horizon, (dS£ m k n / dt)^ w , as 
functions of n. In Figs. EHHl we show the spectrum of ( dEi m kl dt) GW as functions of k. We 
find that the shape of the energy spectrum at the horizon is very similar to that at infinity. 
We thus follow the above procedure to determine the range of summation of k and n for the 
energy spectrum at the horizon. 

For the m-mode, we compute all m-modes from I to —I. For £-mode, we set the maximum 
of £ to be £ max = 20. This value is chosen to reduce computation time. Since the energy flux 
for each £-mode monotonically decreases with increasing £, the relative error of the total flux 
due to this truncation is estimated as 



4.2. Rates of change of the three constants of motion 

We first show the results of the rate of change of energy and the energy spectrum in the 
case of a Schwarzschild black hole. Indeed, the Schwarzschild cases are useful for verifying 
the accuracy of our code. In these cases, the total energy flux after summing the £-, m-, k- 
and n-modes is independent of the inclination angle. We verify this property by computing 
the rate of change of energy when the particle moves on the plane with an inclination angle 
from the equatorial plane around a Schwarzschild black hole. In Table HJ we show the data 
for different inclination angles and eccentricities. We find that, even if the inclination angle 
is changed, (dS / dt)°° remains within the accuracy of 10~ 8 — 10~ 10 . In this table, truncation 




(4-9) 
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Fig. 3. Modal energy flux (dSe m k n I '<^)gw a * infinity for k = and (p, ^inc) — (10M, 20°) in the 
Schwarzschild case. The eccentricities are 0.1 (top left), 0.5 (top right), 0.7 (bottom left) and 
0.9 (bottom right). 



errors A?? <. are shown in square brackets. A?? , is similar to the truncation errors of n- 

('max ) \ m ax / 

and /c-modes, 10~ 10 . Thus, the errors in Table [I] are consistent with the truncation errors of 
the summation of the £-, n- and /c-mode. 

In Table [Til we show similar results for the horizon flux. In these cases, the truncation 
error from the £-mode summation is negligible. The relative error in Table [TTJ is due to the 
truncation of the n- and fc-mode summation. 

In Table IIII| we show the rates of change of the three constants of motion due to the 
emission of gravitational waves to infinity in the case of various bound orbits around a Kerr 
black hole with a = 0.9M. In this table, we use the same orbital parameters as those used 
by Drasco and Hughes in Ref. 21). Our results are consistent with theirs except for the rates 
of change of the Carter constant, (dC/dt). This is because they used formulas for (dC/dt) 
under the approximation that the inclination angle does not change. Of course, this holds 
only approximately. This is the first time that the rate of change of the Carter constant has 
been computed accurately using an adiabatic approximation. In Table. IIIIl we also show 
results for the highly eccentricity of e = 0.9. In Fig. [151 we plot a highly eccentric orbit with 
e = 0.9, p = 6M and 9 mc = 20°. As indicated in the bottom right figure in Fig. HI the peak 
location of the n-mode is approximately 600 when I = m = 20 and k = 0. 

The values in the square brackets in Table UTTI are the truncation errors A?? s, which 
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Fig. 4. Modal energy flux (d£i m k n I <^)gw a * infinity in the Kerr case (a = 0.9M) for k = and 
(p,0 inc ) = (6M,20°). The eccentricities are 0.1 (top left), 0.5 (top right), 0.7 (bottom left) and 
0.9 (bottom right). 



take values from 10~ n to 10 -6 . For highly eccentric cases, e > 0.5, ^(f max ) is larger than 
the truncation error of the n- and /c-mode summation. Thus, the accuracy for such cases is 
only limited by the truncation of the £-mode summation. For e < 0.5, the truncation errors 
of the £-mode summation and n- and fc-mode summation are comparable, and both errors 
contribute to the error of the total flux. 

In Table IPV] we show the time-averaged rates of change of the three constants of motion, 
(d£/dt) n , (dC z /dt) n and (dC/dt) K , due to the absorption at the horizon for various bound 
orbits around a Kerr black hole. The orbital parameters used here are again the same as those 
used by Drasco and Hughes in Ref. 21) except for when e = 0.9. Our results are consistent 
with theirs except for the rate of change of the Carter constant. As in the case of Table ITO| 
this is the first time that the rate of change of the Carter constant due to absorption by 
the black hole has been computed accurately using an adiabatic approximation. The cases 
when e = 0.9 are also new results of this work. The numbers in square brackets are the 
truncation errors of the £-mode summation, 4^), with £ max = 20. The truncation error of 
the £-mode summation is much smaller than that of the n- and /c-mode summation, which 
is approximately 10~ 10 . The accuracy of the data in Table HVl is limited by the truncation 
of the n- and fc-mode summation. 

Here, we discuss the sign of (d£/dt) n . From Table HV] we find that the particle loses 
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Fig. 5. Modal energy flux {d£^ m kl '^)gw a ^ infinity in the Schwarzschild case. (p,e,6i nc ) = 
(10M,0.7,20°). I = 2 (top left), £ = 3 (top right), £ = 4 (bottom left) and I = 5 (bottom 
right). These figures show that the peaks are located at approximately k = 



m. 



energy, i.e., (dS/dt) < 0, when e = 0.3, 0.5 and 0.7 and 9 inc = 80°, whereas the particle gains 
energy, i.e., (dS/dt) n > 0, in the other cases when the eccentricity or the inclination angle is 
small. From Eq. (I2-35H . we find that the sign of each mode {dS/dt)f mkn is determined by the 
sign of aimkn, he., P = oj m kn — ma/ (2Mr + ) = kQo + nfl r -m(a/(2Mr + ) — Q^). The sign 
of (dS/dt) K is determined by the sign of each ai m k n and the absolute value of (d£/dt)f mkn . 
For a > (corotation of the particle and the black hole), when a is large and the modes 
with m > and small k and n dominate the total energy flux, the particle can gain energy. 
We find from Figs. [10] and [13] that when the eccentricity and inclination angle are small, the 
mode with i = m = 2 and k — dominates the total energy flux and the particle gains 
energy. On the other hand, when the eccentricity and inclination angle are large, we find 
from Figs. [TO] and [l~4l that modes with large k and n, which result in P > 0, contribute to 
the total energy flux. Furthermore, Fig. [14] shows that the peak value of m < modes is 
very similar to that for m > 0. This also contributes to making (d£/dt) R positive in the 
cases of large eccentricity and a large inclination angle. 

To confirm the accuracy of the numerical code, we compare our results with the analytical 
post-Newtonian formulas for orbits that are slightly eccentric but highly inclined. 20 '' We show 
the results for p = 100M and various e and 6 inc in Table [V] We find that our numerical 
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Fig. 6. Modal energy flux (d£g m k/dt)Qyy at infinity in the Schwarzschild case, (p, e, #i nc ) = 
(10M,0.7,70°). £ = 2 (top left), £ = 4 (top right), £ = 6 (bottom left) and £ = 8 (bottom 
right). These figures show that the peaks are located at approximately k = £ — m. 



results and the post-Newtonian formulas agree with an accuracy of ~ 10~ 4 or better. Note 
that when p is smaller than 100M, the accuracy of the post-Newtonian formulas becomes 
worse than this value. 

Once we have the rates of change of the constants of motion P = (S, C z , C), we can derive 
the rate of change of the orbital elements i 1 = (p, e, #i nc ). We have the following relation, 

where = dP jdiK In Table [VT| we compare {di l /dt)°° with values derived using the post- 
Newtonian formulas by applying Eq. (14-101) to the data in Table |VJ The relative errors of 
(dp/dt)°° and (de/dt) 00 are 10 -4 , whereas the relative error of {d8 mc /dt)°° is approximately 
10~ 2 . This is because in Ref. 20), (d9 inc / dt)°° is derived only up to 1PN order from the 
leading order, whereas (dp/dt) and (de/dt) are derived up to 2.5PN order. 

We apply Eq. (14-lOp to the data in Tables UTTl and UVI to obtain the rates of change of the 
orbital elements due to the emission of gravitational waves to infinity and to absorption by 
the black hole, which are shown in Tables [VTT1 and IVIIIl respectively. We find that in most 
cases, the flux at infinity and absorption by the black hole exhibit opposite effects except 
when e = 0.3, 0.5 and 0.7 and 9 inc = 80°. In all cases, since the flux at infinity dominates 
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(6M, 0.7, 20°). £ = 2 (top left), I = 3 (top right), I = 4 (bottom left) and £ = 5 (bottom right). 
These figures show that the peaks are located at approximately k = t — m. 

the sign of the total rate of change, the total rates of change shown in Table [IX] have the 
same sign as those for infinity. 

To demonstrate some aspects of the evolution of the orbital elements, in Figs. d6] and [T7] 
we plot the total rates of change ((dp/dt) , (d9 inc /dt)) and ((dp/dt) , (de/dt)) on the (p, # inc ) 
and (p, e) planes, respectively, in the case when a = 0.9M and £ max = 5. We find that 
although the eccentricity is always decreasing at large p, it can increase near the last stable 
orbit (LSO). The change in the inclination angle is not very significant in the figures, but it 
is always increasing at large p. 
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8. Modal energy flux {d8^ m k/ 'dt)^^ at infinity in the Kerr case (a = 0.9M). (p, e, 0i nc ) = 
(6M, 0.7, 80°). 1 = 2 (top left), I = 4 (top right), £ = 6 (bottom left) and £ = 8 (bottom right). 
These figures show that the peaks are located at approximately k = £ — m. 
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9. Modal energy flux | (df/ m fc„/(It) GW | at the horizon for k = and (p, 9i nc ) = (10M, 20°) in 
the Schwarzschild case. The eccentricities are 0.1 (top left), 0.5 (top right), 0.7 (bottom left) 
and 0.9 (bottom right). 
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Fig. 10. Modal energy flux | (d£g m k n /dt}^,^ | at the horizon in the Kerr case (a = 0.9M) for k = 
and (p,e inc ) = (6M,20°). The eccentricities are 0.1 (top left), 0.5 (top right), 0.7 (bottom left) 
and 0.9 (bottom right). 
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Fig. 11. Modal energy flux | {dSg m k/ dt)^^ | at the horizon in the Schwarzschild case. (p,e,0i nc ) = 
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Fig. 12. Modal energy flux | (dS^^/dt)^^ | at the horizon in the Schwarzschild case, (p, e, 6[ nc ) = 
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Fig. 13. Modal energy flux 
(p,e,0 inc ) = (6M,0.7,20°) 
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(d£e m k/dt) G yj | at the horizon in the Kerr case (a = 0.9M). 
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(bottom right). These figures show that the peaks are located at approximately k 
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Fig. 14. Modal energy flux | (d£g m k / 'dt) Gw | at the horizon in the Kerr case (a = 0.9M). 
(p,e,0 inc ) = (6M,0.7,80°). £ = 2 (top left), I = 4 (top right), £ = 6 (bottom left) and I = 8 
(bottom right). These figures show that the peaks are located at approximately k = 



m. 




Y X' 

Fig. 15. Plots of the orbits in the same coordinate systems as those used in Fig. [2j This generic 
geodesic orbit has eccentricity e = 0.9, semilatus rectum p = 6M and inclination angle 6 mc = 
20°. The spin of the black hole is set to a = 0.9M. 
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Table I. Time- averaged rates of change of the energy of a particle due to the emission of gravi- 
tational waves to infinity for the case of a Schwarzschild black hole. In this table, the orbital 
radius is 10M. We compare the result for the equatorial plane with that for a nonequatorial 
plane. Truncation errors A?? A are shown in square brackets. Here we set ^ max = 20. These 
results show that the error of {d£ / dt)°° is consistent with the truncation error of the mode 
summation. 



a/M 


p/M 


e 


^inc 


(de/dt)°° 


Relative error 





ID 


0.1 


no 




-6.31752474720 x 10' 


-5 


M n — 111 

[10 iX J 





10 


0.1 


20° 


-6.31752474730 x 10" 


-5 


1.6xl0- n 





10 


0.1 


45° 


-6.31752474742 x 10' 


-5 


3.5xl0- n 





10 


0.1 


70° 


-6.31752474665 x 10" 


-5 


8.7xl0- n 





10 


0.5 


0° 


-9.27335011503 x 10" 


-5 


[1CT 11 ] 





10 


0.5 


20° 


-9.27335011442 x 10" 


-5 


6.6xl0- n 





10 


0.5 


45° 


-9.27335011373 x 10" 


-5 


1.4X10" 10 





10 


0.5 


70° 


-9.27335011191 x 10" 


-5 


3.4xl0~ 10 





10 


0.7 


0° 


-9.46979134409 x 10" 


-5 


[io- 9 ] 





10 


0.7 


20° 


-9.46979134028 x 10" 


-5 


4.0xl0~ 10 





10 


0.7 


45° 


-9.46979133931 x 10" 


-5 


5.0xl0" 10 





10 


0.7 


70° 


-9.46979131018 x 10" 


-5 


3.6xl0~ 9 





10 


0.9 


0° 


-4.19426469206 x 10" 


-5 


[io- 8 ] 





10 


0.9 


20° 


-4.19426468442 x 10" 


-5 


1.8xl0~ 9 





10 


0.9 


45° 


-4.19426468407 x 10" 


-5 


1.9xl0- 9 





10 


0.9 


70° 


-4.19426437158 x 10" 


-5 


7.6xl0~ 8 
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Table II. Time- averaged rates of change of the energy of a particle due to the emission of gravita- 
tional waves to the horizon for the case of a Schwarzschild black hole. In this table, the orbital 
radius is 10M. We compare the result for the equatorial plane with that for a nonequatorial 
plane. Truncation errors A,„ n are shown in square brackets. Here we set £ max = 20. 

(.{-max ) 



a/M p/M e 6 mc 


(dS/dtf 


Relative error 


10 0.1 0° 
10 0.1 20° 
10 0.1 45° 
10 0.1 70° 


-1.53365819445 x IO" 8 
-1.53365819444 x IO" 8 
-1.53365819443 x 10~ 8 
-1.53365819444 x IO" 8 


[10~ 35 ] 
6.5xl0~ 12 
1.3xl0- n 
6.5xl0~ 12 


10 0.5 0° 
10 0.5 20° 
10 0.5 45° 
10 0.5 70° 


-1.41298859260 x 10~ 7 
-1.41298859246 x 10~ 7 
-1.41298858862 x 10~ 7 
-1.41298859240 x 10^ 7 


[io- 31 ] 

9.9xl0~ n 
2.8xl0~ 9 
1.4xl0^ 10 


10 0.7 0° 
10 0.7 20° 
10 0.7 45° 
10 0.7 70° 


-3.55415030114 x 10~ 7 
-3.55415029914 x IO" 7 
-3.55415027484 x 10~ 7 
-3.55415029916 x 10~ 7 


[io- 26 ] 

5.6xl0- 10 
7.4xl0~ 9 
5.6xl0" 10 


10 0.9 0° 
10 0.9 20° 
10 0.9 45° 
10 0.9 70° 


-3.65214284306 x IO" 7 
-3.65214284171 x 10~ 7 
-3.65214284814 x 10~ 7 
-3.65214285085 x 10~ 7 


[io- 22 ] 

3.7xl0~ 10 
1.4xl0~ 9 
2.1xl0- 9 
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Table III. Time-averaged rates of change of the three constants of motion, energy {d£/dt)°° , 
angular momentum {dC z /dt)°° and the Carter constant {dC/dt)°°, due to gravitational waves 
radiated to infinity per unit mass for various generic orbits around a Kerr black hole. The 
orbital parameters used here are the same as those used by Drasco and Hughes in Ref. 21) 
except when e = 0.9. Our results are consistent with theirs except for the rate of change of 
the Carter constant, (dC/dt)°°. Numbers in square brackets are the truncation errors of the 
■£-mode summation, Z\(^ max ) with £ mSlX = 20. Note that the case of a = 0.9M, p = 6M, e = 0.9 
and 6 mc = 80° does not result in stable bound orbits. 



a/M 


p/M 


e 


^inc 


(dS/dt)°° 






(dC/dt)°° 




0.9 


6 


0.1 


20° 


-5.87363800087 x 10- 4 [10- u ] 


-8.53727881580 x 10~ 


-3 


-5.24019848546 x 10" 


-3 


0.9 


6 


0.1 


40° 


-6.18322941497 x lO^lO" 11 ] 


-7.63099401313 x 10~ 


-3 


-2.02271137874 x 10" 


-■> 


0.9 


6 


0.1 


60° 


-6.83348195277 x lO^lO" 11 ] 


-6.07829111698 x 10~ 


-3 


-4.32194650210 x 10" 


-■> 


0.9 


6 


0.1 


80° 


-8.05858117692 x 10- 4 [10- 10 ] 


-3.62538058308 x 10" 


-3 


-7.18520476701 x 10" 


-■> 


0.9 


6 


0.3 


20° 


-6.80409929713 x 10- 4 [10- 9 ] 


-8.62590762129 x 10" 


-3 


-5.22145052503 x 10" 


-3 


0.9 


6 


0.3 


40° 


-7.26541780924 x lO^lO" 9 ] 


-7.84019287653 x 10" 


-3 


-2.04389407707 x 10" 


-2 


0.9 


6 


0.3 


60° 


-8.30597576584 x 10- 4 [10- 9 ] 


-6.49674204013 x 10" 


-3 


-4.50701803710 x 10" 


-2 


0.9 


6 


0.3 


80° 


-1.08394107072 x 10- 3 [HT 9 ] 


-4.38279817141 x 10" 


-3 


-8.18315782169 x 10' 


-2 


0.9 


6 


0.5 


20° 


-7.98925629079 x 10- 4 [10- 8 ] 


-8.34750401557 x 10" 


-3 


-4.94704500000 x 10" 


-3 


0.9 


6 


0.5 


40° 


-8.74335722008 x 10- 4 [10- 8 ] 


-7.81941309824 x 10" 


-3 


-1.98900691857 x 10' 


-■> 


0.9 


6 


0.5 


60° 


-1.05884649558 x 10- 3 [10- 8 ] 


-6.95065613502 x 10" 


-3 


-4.65662839707 x 10" 


-■> 


0.9 


6 


0.5 


80° 


-1.67699406035 x 10- 3 [10- 7 ] 


-5.90867286937 x 10" 


-3 


-1.01809636563 x 10" 


-1 


0.9 


6 


0.7 


20° 


-7.73126177805 x lO^lO" 7 ] 


-6.69310061924 x 10" 


-3 


-3.88688412163 x 10" 


-3 


0.9 


6 


0.7 


40° 


-8.75195550414 x 10- 4 [10- 7 ] 


-6.53053154749 x 10" 


-3 


-1.62422254375 x 10" 


-2 


0.9 


6 


0.7 


60° 


-1.14691287812 x lO^lO" 7 ] 


-6.38027895400 x 10" 


-3 


-4.14632893132 x 10" 


-2 


0.9 


6 


0.7 


80° 


-2.71933022663 x lO^lO" 6 ] 


-8.40183718587 x 10" 


-3 


-1.34155501090 x 10" 


-1 


0.9 


6 


0.9 


20° 


-3.22243369277 x 10- 4 [10- 6 ] 


-2.36914164090 x 10" 


-3 


-1.35698436703 x 10" 


-3 


0.9 


6 


0.9 


40° 


-3.81407518944 x 10- 4 [10- 6 ] 


-2.43273577381 x 10" 


-3 


-5.96550294295 x 10" 


-3 


0.9 


6 


0.9 


60° 


-5.55433712227 x lO^lO" 6 ] 


-2.68099758625 x 10" 


-3 


-1.71022011422 x 10" 


-2 


0.9 


6 


0.9 


80° 
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Table IV. Time-averaged rates of change of the three constants of motion, energy (d£/dt) , 
angular momentum (dC z /dt) and the Carter constant (dC/dt}^ 1 due to gravitational waves 
absorbed at the horizon per unit mass for various generic orbits around a Kerr black hole. The 
orbital parameters used here are the same as those used by Drasco and Hughes in Ref. 21) 
except when e = 0.9. Our results are consistent with theirs except for the rate of change of the 
Carter constant, (dC/dt)^. Numbers in square brackets are the truncation errors of the ^-mode 
summation, Z\(£ max ) with ^ max = 20. The case of q = 0.9M, p = 6M, e = 0.9 and 8\ nc = 80° 
does not result in stable bound orbits. 



a/M 


p/M 


e 


#inc 


(d£/dt) H 


(d£ z /dt) H 




(dC/dt) H 




0.9 


6 


0.1 


20° 


4.25245612585 x 10- 6 [10- 25 ] 


6.71500254679 x 10" 


-5 


1.41073696418 x 10" 


6 


0.9 


6 


0.1 


40° 


3.94882721384 x lO^lO" 27 ] 


7.73637887857 x 10" 


-5 


-2.22131838862 x 10 


-5 


0.9 


G 


0.1 


60° 


3.33113477148 x 10- 6 [10- 30 ] 


1.11677030642 x 10" 


-4 


-1.02353953718 x 10 


-4 


0.9 


6 


0.1 


80° 


9.50601680011 x 10- 7 [10- 35 ] 


1.90137812316 x 10" 


-4 


-2.72562414552 x 10 


-4 


0.9 


6 


0.3 


20° 


5.86967815445 x lO^lO" 23 ] 


7.76727457985 x 10" 


-5 


-5.86205174377 x 10 


-6 


0.9 


G 


0.3 


40° 


5.84180692574 x lO^lO" 23 ] 


1.00066438717 x 10" 


-4 


-5.95214221681 x 10 


-5 


0.9 


G 


0.3 


60° 


5.19494530315 x 10- 6 [10- 26 ] 


1.66116835805 x 10" 


-4 


-2.26658646533 x 10 


-4 


0.9 


6 


0.3 


80° 


-2.95984810527 x 10- 9 [10- 22 ] 


3.45153570099 x 10" 


-4 


-7.00228568783 x 10 


-4 


0.9 


6 


0.5 


20° 


8.34425799664 x lO^lO" 21 ] 


9.13622955381 x 10" 


-5 


-2.13883155844 x 10 


-5 


0.9 


G 


0.5 


40° 


8.94532600748 x 10- fi [10- 21 ] 


1.37186174930 x 10~ 


-4 


-1.42441848989 x 10 


-4 


0.9 


6 


0.5 


60° 


8.08290137218 x lO^lO" 23 ] 


2.69931271416 x 10" 


-4 


-5.34259432172 x 10 


-4 


0.9 


6 


0.5 


80° 


-5.98308999615 x 10- 6 [10- 17 ] 


7.40794676040 x 10" 


-4 


-2.07799385090 x 10" 


-3 


0.9 


6 


0.7 


20° 


9.29526284834 x 10- 6 [10- 19 ] 


8.90473867250 x 10~ 


-5 


-3.96683923123 x 10 


-5 


0.9 


6 


0.7 


40° 


1.05570527008 x 10- 5 [10- 19 ] 


1.56425173930 x 10~ 


-4 


-2.43302981071 x 10 


-4 


0.9 


G 


0.7 


60° 


9.20041590887 x lO^lO" 21 ] 


3.62722481332 x 10" 


-4 


-9.65381244057 x 10 


-4 


0.9 


6 


0.7 


80° 


-2.74019070298 x 10- 5 [10- 13 ] 


1.59338857620 x 10" 


-3 


-5.85248852167 x 10~ 


-3 


0.9 


6 


0.9 


20° 


4.17773280477 x 10- 6 [10- 18 ] 


3.74778275669 x 10" 


-5 


-2.72898298129 x 10 


-5 


0.9 


G 


0.9 


40° 


4.94483767585 x lO^lO" 18 ] 


7.68457866865 x 10" 


-5 


-1.65632369297 x 10 


-4 


0.9 


G 


0.9 


60° 


3.94007307443 x 10- 6 [10- 14 ] 


2.10536770888 x 10" 


-4 


-7.17182301393 x 10 


-4 


0.9 


6 


0.9 


80° 
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Table V. Comparison of the time-averaged rates of change of the three constants of motion derived 
using our numerical method and analytical post-Newtonian expressions 20 -* for orbits that are 
slightly eccentric but greatly inclined in the case of a = 0.9M and p = WOM. Our numerical 
results are consistent with the post-Newtonian results. Relative errors are always approximately 
10" 4 . 



e 


^inc 


(df/*) Numerical 


(dCJdt) 


OO 

Numerical 


(dC/ ^Numerical 




—Newton 


{ dC z / dt) Post _ Newton 


(dC/di)~ st _ Newton 


0.01 


20° 


-6.2072 x 10" 


10 


-5.8382 x 10" 


-7 


-1.4687 x 10" 


-0 


-6.2067 x 


10" 


10 


-5.8376 x 10" 


-7 


-1.4687 x 10" 


-6 


0.01 


45° 


-6.2142 x 10" 


10 


-4.4029 x 10" 


-7 


-6.2905 x 10" 


-6 


-6.2136 x 


10" 


-10 


-4.4024 x 10" 


-7 


-6.2901 x 10" 


-6 


0.01 


70° 


-6.2259 x 10 


-10 


-2.1420 x 10" 


-7 


-1.1146 x 


10" 




-6.2251 x 


10" 


-10 


-2.1415 x 10" 


-7 


-1.1145 x 10" 


-5 


0.05 


20° 


-6.2302 x 10" 


10 


-5.8299 


x 10" 


-7 


-1.4667 x 


10" 


-6 


-6.2296 x 


10" 


10 


-5.8293 x 10" 


-7 


-1.4666 x 10" 


-6 


0.05 


45° 


-6.2373 x 10" 


10 


-4.3967 


x 10" 


-7 


-6.2816 x 


10" 


-6 


-6.2366 x 


10" 


-10 


-4.3962 x 10" 


-7 


-6.2812 x 10" 


-6 


0.05 


70° 


-6.2490 x 10" 


10 


-2.1390 


x 10" 


-7 


-1.1130 x 


10" 




-6.2481 x 


10" 


-10 


-2.1386 x 10" 


-7 


-1.1129 x 10" 


-5 


0.09 


20° 


-6.2827 x 10" 


10 


-5.8103 


x 10" 


-7 


-1.4617 x 


10" 


-6 


-6.2818 x 


10" 


10 


-5.8097 x 10" 




-1.4616 x 10" 


-6 


0.09 


45° 


-6.2899 x 10" 


10 


-4.3820 x 10" 


-7 


-6.2606 x 


10" 


-6 


-6.2890 x 


10" 


-10 


-4.3814 x 10" 




-6.2601 x 10" 


-6 


0.09 


70° 


-6.3019 x 10 


-10 


-2.1320 


x 10" 


-7 


-1.1093 x 


10" 




-6.3008 x 


10" 


10 


-2.1315 x 10" 


-7 


-1.1092 x 10" 


-5 



Table VI. Comparison of the time-averaged rates of change of orbital elements derived using our 
numerical method and the analytical post-Newtonian expressions 20 ^ for orbits that are slightly 
eccentric but greatly inclined in the case of a = 0.9M and p = lOOAf . Our numerical results 
are consistent with the post-Newtonian results. The relative errors of both (dp/dt) and (de/dt) 
are always approximately 10 -4 . However, the relative errors of (dQ- mc /dt) are approximately 

lO" 2 Since <^incM>Post-Newton in Ref ' 20 ) is 1PN - 



e 




{dp/ ^Numerical 


(de/dt) ^^1 


{d0 iuc /dt) 


OO 

Numerical 


(dp/dt)^ 


-Newton 


(de/ dt)p ost 


-Newton 


{ d9i nc / dt) Post-Newton 


0.01 


20° 


-1.2560 x 10~ 5 


-1.9807 x 


10" 




2.5353 x 10 


-9 


-1.2558 x 10" 


-5 


-1.9803 x 


10 


-9 


2.6807 x 10~ 9 


0.01 


45° 


-1.2587 x lO" 5 


-1.9851 x 


10" 


-9 


5.3432 x 10 


-9 


-1.2585 x 


10" 




-1.9846 x 


10 


-9 


5.5488 x lO" 9 


0.01 


70° 


-1.2631 x Hr 5 


-1.9923 x 


10" 


-9 


7.3134 > 


10 


-9 


-1.2630 x 


10" 


-5 


-1.9917 x 


10 


-9 


7.3877 x lO" 9 


0.05 


20° 


-1.2541 x lO" 5 


-9.8778 x 


10" 




2.5450 > 


io- 


-9 


-1.2540 x 


10" 


-5 


-9.8753 x 


10 


-9 


2.6900 x 10" 9 


0.05 


45° 


-1.2568 x 10" 5 


-9.8996 x 


10" 


-9 


5.3636 > 


io- 


-9 


-1.2567 x 


10" 


-5 


-9.8968 x 


10" 


-9 


5.5683 x 10~ 9 


0.05 


70° 


-1.2613 x Hr 5 


-9.9357 x 


10" 


-9 


7.3412 > 


io- 


-9 


-1.2611 x 10" 


-5 


-9.9324 x 


10 


-9 


7.4142 x 10~ 9 


0.09 


20° 


-1.2497 x lO" 5 


-1.7671 x 10" 


-8 


2.5672 x 10- 


-9 


-1.2496 x 10" 


-5 


-1.7666 x 10" 


-8 


2.7112 x 10" 9 


0.09 


45° 


-1.2524 x lO" 5 


-1.7711 x 


10" 


-8 


5.4104 > 


io- 


-9 


-1.2523 x 


10" 


-5 


-1.7705 x 


10- 


-8 


5.6127 x 10~ 9 


0.09 


70° 


-1.2569 x 10- 5 


-1.7775 x 10" 


-8 


7.4050 > 


10 


-9 


-1.2567 x 10" 


-5 


-1.7768 x 


io- 


-8 


7.4748 x 10" 9 



§5. Summary 

In this paper, we developed a numerical code to compute gravitational waves induced 
by a particle orbiting around a Kerr black hole. We obtained the rates of change of energy, 
angular momentum and the Carter constant for various eccentric and inclined orbits. This 
is the first time that the rate of change of the Carter constant has been evaluated accurately. 
These computations include highly eccentric cases, i.e., e = 0.9. In previous works, such 
high eccentricity was not treated. 

Our numerical method is mainly divided into four parts: the computation of the radial 
and polar motion, the homogeneous solution of the Teukolsky equation, the integration of 
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Fig. 16. Evolution of eccentric and inclined orbits projected on the e = const, plane. The black 
hole spin is a = 0.9M. The dashed curves represent the last stable orbit (LSO). Each point rep- 
resents an orbit and each arrow represents the rate of change of the orbit ((dp/dt) , (d6- mc /dt)). 
Here, the lengths of arrows are normalized appropriately. The top left and top right figures 
show the evolution on the e = 0.1 and e = 0.3 planes respectively. The bottom left and bottom 
right figures show the evolution on the e = 0.5 and e = 0.7 planes respectively. At a large 
distance, p is always decreasing and 9- mc is always increasing. In the computation of this figure, 
we set £ max = 5. 



Eq. f l2-32p and the mode summation in Eqs. (I2-35H ( I2-38H . We found that the radial and 
polar motion can be expressed in terms of Jacobi elliptic functions. This enabled us to solve 
the geodesic motion more accurately than the method in which the equations of motion 
are integrated numerically. The asymptotic amplitudes, such as B^^, can be computed 
directly from the MST formalism. We do not need to evaluate the homogeneous solutions at 
a very large distance to obtain the asymptotic amplitudes. We computed the homogeneous 
Teukolsky solution, R^ muJ , at a radius r between r min < r < r max using the MST formalism 
accurately. The homogeneous solution at other radial points is computed by successive Taylor 
series expansions. The Taylor series method gives very high accuracy and is much faster than 
the use of the hypergeometric function expansion at all radial points. We computed the 
integral, Eq. (12-321) . using the trapezium rule, which gives very high accuracy when periodic 
functions are integrated over one period. 

The accuracy of the numerical results are limited by the truncation of the summation of 
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Fig. 17. Evolution of eccentric and inclined orbits projected on the 6 inc = const, plane. The black 
hole spin is a = 0.9M. The dashed curves represent the last stable orbit (LSO). Each orbit 
is represented by a point and its evolution is represented by a vector ((dp/dt) , (de/dt)). The 
top left figure shows the evolution on the 9- inc = 20° plane and the top right figure shows the 
evolution on the inc = 40° plane. The bottom left figure shows the evolution on the 9- mc = 60° 
plane and the bottom right figure shows the evolution on the 6 inc = 80° plane. At a large 
distance, p and e are always decreasing. In the computation of this figure, we set £ max = 5. 



the k- and n-modes. We have verified the behavior of the energy spectrum of the k- and 
n-modes. We determined the range of the summation of k and n to obtain an error due to 
the truncation of k- and n-modes of less than 10 -10 . We truncated the £-mode at I = 20. 
This value was chosen to reduce computation time. The error due to the truncation of the 
^-mode depends on the orbital parameters. When the eccentricity is small, i.e. e < 0.3, this 
error is approximately 10 -9 . However, when the eccentricity is e = 0.9, this error becomes 
10~ 5 , which is the largest error out of the results computed in this paper. Note that since 
the error is only limited by the truncation of the £-, k- and n-modes, it is straightforward to 
improve the accuracy. 

To confirm the accuracy of our code, we computed the energy flux from a Schwarzschild 
black hole for cases when the orbits are inclined with respect to the equatorial plane. In 
the Schwarzschild case, the energy flux should not depend on the inclination angle. Thus, 
we can estimate the accuracy of the code by comparing the results for inclined orbits and 
equatorial orbits. We found that the accuracy of our code is consistent with the estimates 
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Table VII. Time- averaged rates of change of orbital elements, semilatus rectum (dp/dt) 00 , ec- 
centricity (de/dt)°° and inclination angle {d9i nc /dt)°° due to gravitational waves radiated to 
infinity per unit mass for various generic orbits around a Kerr black hole. The orbital param- 
eters used here are the same as those used in Tables Hill and HVl Here we set £ max = 20. Note 
that the case of q = 0.9M, p = 6M, e = 0.9 and 9 inc = 80° does not result in stable bound 
orbits. 
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of the truncation errors of the k- and n-modes. 

Although we have not completely optimized the code, we briefly discuss the computation 
time here. In the current code, the computation time for one mode is roughly 0.1 — 0.3 
seconds. The total computation time is about 6 — 12 hours when q = 0.9, p = 6M, e = 0.1 
and 6 mc = 20° - 80°, about 1 day when q = 0.9, p = 6M, e = 0.7 and 6 mc = 20° - 80°, 
and 1 — 3 days when q = 0.9, p = 6M,e = 0.9 and 9 inc = 20° — 60°. These times are for 
the results of computation with one 2.5 GHz AMD Opteron CPU. In Fig.l in Ref. 14), the 
computation time for one mode is shown. The computation time using 8 CPUs is about 0.4 
seconds when q = 0.7, p = 10M, e = 0.5 and 6- nic = 45°. Thus, the computation time using 
one CPU will be a few seconds. Thus, our computation time appears to be much shorter 
than that in Ref. 14). 

By optimizing the code, we can further increase its speed. Also, it is possible to consider 
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Table VIII. Time-averaged rates of change of orbital elements, semilatus rectum (dp/dt) , eccen- 
tricity (de/dt) and inclination angle (d8- mc /dt)^ due to absorption by a black hole per unit 
mass for various generic orbits around a Kerr black hole. The orbital parameters used here are 
the same as those used in Tables IIIII and IIVI Here we set l mSuX = 20. Note that the case of 
q = 0.9M, p = 6M, e = 0.9 and 9 inc = 80° does not result in stable bound orbits, (dp/dt)^ 
and (dOinc/dt) have opposite signs to (dp/dt)°° and (d6\ nc / dt)°° respectively, (de/dt) also 
has the opposite sign to (de/dt)°° except when 9- mc = 80°. 
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variants of the method used to obtain the homogeneous solutions, the orbital motion, and so 
forth. Although we have not used numerical integration methods, it may be advantageous 
to use them under some circumstances. We will consider these issues in the future. We will 
also investigate the computation of gravitational waves including the effects of the adiabatic 
evolution of a particle orbit due to the emission of gravitational waves. This is important for 
investigating EMRI through the data analysis of LISA. We will also investigate the possibility 
of computation over a wider range of orbital parameters, such as an eccentricity larger than 
0.9, and an inclination angle of 9 = 90°. 

We are planning to make our code available to the public so that a wide range of people 
may use it. We believe our code will be very useful for investigating astrophysical and data 
analysis issues of EMRI and for analyzing the data obtained from LISA, DECIGO and BBO. 
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Table IX. Total time-averaged rates of change of orbital elements, semilatus rectum (dp/dt), 
eccentricity (de/dt) and inclination angle {d9{ nc /dt) due to gravitational waves per unit mass 
for various generic orbits around a Kerr black hole. The values in this table are the sums of 



the values in Tables IVIII and IVIIII 
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Appendix A 

Explicit Expressions for Z^ nuJ and Zf^ 

In this section, we give explicit expressions for the amplitude of the partial waves zfj^ '. 
Using the Green function of the radial Teukolsky equation, the solutions of the Teukolsky 
equation are expressed as 

RemUr) = + R? m Jr)Z? rnuj ( r ), (A-l) 
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where R™J^ '(r) satisfy ingoing/outgoing wave conditions at the horizon/infinity, and 
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Appendix B 

Geodesic Motion in the Kerr Spacetime 



In this section, we discuss the solutions of the geodesic equations in detail. First, we 
describe analytical solutions of the r and 9 components of the geodesic equations, which are 
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expressed as 



= (Bl) 
= 0{cos6), (B-2) 



where 



P(r) = £(r 2 + a 2 ) - aC z , 
R(r) = [P(r)] 2 - A[r 2 + (a£ - C z f + C], 
0(cos9) = C-{C + a 2 (l - £ 2 ) + £ 2 ) cos 2 9 + a 2 (l - E 2 ) cos 4 0. 

Since i?(r) and <9(cos#) are fourth-order polynomials, there are four zeros of r and cos# 
for each function, respectively. A geodesic can be specified if we set two zero points, r min and 
r max; for the radial part and one zero point, cos# mm , for the polar part. This corresponds to 
the fact that a one-to-one correspondence exists between (r min , r max , 6* min ) and (E,C Z ,C). 

It is convenient to introduce the orbital parameters, eccentricity e, semilatus rectum p 
and inclination angle 9- mc , defined as 

'"min Z \ i '"max Z > ^inc ~\~ (sgn/^) $min ~Z- (B '3) 

1 + e 1 — e 2 

The three constants of motion, (E, C z , C), are expressed in terms of these orbital parameters 
(p,e,9 mc ) 2 ^ 

To solve the differential equations for r and cos^, we rewrite R(r) and 6>(cos6>) as 

R(r) = (1 — £ 2 )(ri — r)(r — T2}{r — t$){t — r$), 
(cos 9) = C 2 z e (z^ — cos 2 9)(z + — cos 2 9), 

where 



p p (A + B) + J (A + B) 2 - AAB AB 
n = —' r2 = —e r3 = 2 ' r4 = 77 ' 

A + B-^-fr + r*), AB= JT ^L_, (B-4) 

and e = a 2 (l - E 2 )/C 2 Z , z_ = cos 2 9 min , z + = C/(C 2 z e z_). 

Let the solutions of Eqs. ( 1B-1I) and ( 1B-2I) in terms of r or be X^ r '(r) and A^(0), 
respectively. The functions \^(r) and A^(#) are expressed as 

2Aq ; (ri) - A^(r) r : n -> r 2 , 
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A (8) 8 : | — > mill , 



A^(0) = <{ 2A^(0 mill )-Aj» 8 
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and the function F is an elliptic integral of the first kind. In the following, we describe the 
elliptic integrals and functions using the notation in Ref. 30). 

The orbital frequencies of the radial and polar motion with respect to A, which are 
denoted by T r and Tg, respectively, are defined in (12-181) . They are expressed as 



r _ *y(l ~£ 2 )(ri ~r 3 )(r 2 - rl) _ 7r£ zv /e^T 

2K(k r ) ' 6 2K(k e ) ' 1 ' 

Here K(k) is the complete elliptic integral of the first kind, and 



V ri - r 3 r 2 - r 4 V z + 

Furthermore, Eq. (1B-5I) can be solved inversely. 

/ \ r 3 {ri -r 2 )sn 2 ((f r (w r );k r ) -r 2 {r 1 - r 3 ) 

= —7 v — 2 , 7 v r v 7 r— , cos^ (w e ) = y/zZm{(pe{w B )\ kg). (B-9) 

(ri - r 2 ) sn^(v9 r (w r ); fc r J - (n - r 3 ) 

For convenience, we have introduced the angle variables, 

W r = T r X, Wg = T e \. (B-10) 
The function sn(</?, k) (and cn(<^, fc) and dn(y9, fc) below) is a Jacobi elliptic function and 



K(k r ) 



(0<W r < 7T) 



(27T- Wr )^, (7T < tlV < 27T) 
(0 < Wg < f) 

^M = <! (tt-^)^, (|<«;,<f) (B-ll) 
(w e -2n) 2 -^l. (f <™ e <7r) 
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By differentiating r and cos# with respect to A, we respectively obtain dr/dX and dcos6/dX 
analytically, which are expressed as 



dr 
dX 



dcos9 



where 



dX 



(w r ) =2sgn 



dip r 
dw. 



sn(ip r ; k r )cn(Lf r ; k r )dn(Lf r ; k r )(r 2 - r 3 )(ri - r 3 )(ri - r 2 ) K(k r )T r 



((ri - r 2 ) sn 2 (v9 r ; k r ) - (n - r 3 )) 2 
=2sgn ( v /iTcn(v9 9 ; /c 9 )dn(y2 ; k e ) K 



dwc 



It 



sgn(xj 



1 for x > 0, 
for x = 0, 
— 1 for x < 0. 
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